1 Install R package, RStudio, directory set-up, GitHub

1.1 Install R and RStudio

Install R and install RStudio https://posit.co/download/rstudio-desktop/ based on your operating system, available for Linux, MacOS, and Windows

1.2 Differential expression analysis (DEA) work directory setting up

Open up RStudio and setting up a new project for DEA analysis 1.Go to the File menu and select New Project.2.In the New Project window, choose New Directory. Then, choose New Project. Name your new directory RNA_seq_DEA and then “Create the project as subdirectory of:” the Desktop or OneDrive - University of Helsinki (or location of your choice).3.The new project should automatically open in RStudio.

1 To check whether or not you are in the correct working directory, use getwd(). The path Desktop/RNA_seq_DEA should be returned to you in the console.

getwd()
## [1] "C:/Users/timoi/OneDrive - University of Helsinki/ZilanWen Github/RNA-seq_Differential-Gene-Expression-Analysis"

To get the location of working directory For example: At Desktop, build a folder At RStudio –> Session –> Set Working Directory –> Choose Directory –> Desktop –> RNA-seq_DEA

or use code below:

setwd("C:/Users/username/Desktop/RNA-seq_DEA")

Check the current R version. The current R version is R-4.4.0

R.Version()#the code in the code chunk will not be run.

Within your working directory use the New folder button in the bottom right panel to create two new directories data and results.

Now we need to download the files from Moodle and the data list which looks like below
data Then we go to the File menu and select New File, then select R Script. This should open up a script editor in the top left hand corner. This is where we will be typing and saving all commands required for this analysis. Now save the file as DEA_script.R. When finished your working directory should now look similar to this:
directories data for the original data, results for generated files and figures.

1.3 GitHub (optional)

GitHub is an online software development platform. It’s used for storing, tracking, and collaborating on software projects. It makes it easy for developers to share code files and collaborate with fellow developers on open-source projects.

  • Set up your GitHub: Go here https://github.com/, log in or sign up with your email.

  • Create your repository: Repository is the project you can create which is on the top left of your GitHub home page.To create a new repository, you can click the green button New, then name the new repository.Choose Public or Private, and initialize this repository with add a README file. Next Click Create repository. The list of repositories is on your Github home page.

  • Use the Github desktop version: We are going to download and upload by commiting changes from our local files to the cloud. Go to GitHub Desktophttps://desktop.github.com/, download based on your operating system. Then go to sign in with your GitHub account.Go back to your GitHUb repository, and click open with Github Desktop
    click open with Github Desktop

  • Clone a repository and connect repository to local files In local path, choose where we want to save the repository in our local computer.click Clone. Now we have cloned the repository that we created on GitHub and we have connected it to the GitHub desktop application.
    1

  • upload local files to GitHub cloud Any changes and updated files can be uploaded to GitHUb cloud though GitHub Desktop.Atfter change files, go back to GitHub Desktop and navigate the corresponding repository, you will have the changes listed on the top left. We now give the change a title such as ‘11’ in thesummary (required)box and click Commit to main on the bottom left, Then on the right, we can push comments to the origin remote by clicking Push origin.
    local path

  • Clone a repository from other GitHUb accounts Forking a repository can make a copy of a repository from other account to your account.

local path
local path

2 Count normalization

Created on 14.05.2024
@author: Zilan Wen
Email:

2.1 Normalization Overview

The first step in the DE analysis workflow is count normalization, which is necessary to make accurate comparisons of gene expression between samples.

The main factors often considered during normalization are:

  • sequencing depth or coverage.
local path
local path

  • Gene length

  • RNA composition

These three factors have been described very well in [DEG workshop Salmon] (https://github.com/hbctraining/DGE_workshop_salmon/blob/master/lessons/02_DGE_count_normalization.md)

2.1.1 Common normalization methods

The differences between RPKM, FPKM and TPM are explained clearly by SatQuest video

2.1.2 DESeq2-normalized counts: Median of ratios method

This method is explained clearly in DEG_workshop_salmon Since tools for differential expression analysis are comparing the counts between sample groups for the same gene, gene length does not need to be accounted for by the tool. However, sequencing depth and RNA composition need to be taken into account.

2.2 Count normalization of our dataset with DESeq2

For transcriptome study, RNA samples which are isolated from control samples and treatment samples are sent to sequencing company. The following steps are :

  • RNA-seq with a sequencing depth of 10-30 M reads per library (at least 3 biological replicates per sample)

  • Aligning or mapping the quality-filtered sequenced reads to respective genome (e.g. HISAT2 or STAR or Bowtie2).

  • Quantifying reads that are mapped to genes or transcripts (e.g.HTseq)

  • Gene differential expression analysis

  • Functional analysis such as enrichment of known biological functions, interactions, or pathways

RNA samples Heterobasidion annosum is a forest pathogen which can cause root and stem rot disease on Scots pine and Norway spruce. In this experiment, two-month old Scots pine plants were inoculated with H. annosum. After one month infection, total RNA was extracted from Scots pine roots. Control plants were fungus-free. We have three biological replicates for each treatment. The treatment samples are Ha1,Ha2, Ha3 and the control samples are Ctr1, Ctr2, Ctr3. Since genome of Scots pine is not available, we align reads to Pinus taeda genome which can be found in plantgenie https://plantgenie.org/.

In this tutorial we start with raw reads count ´.txt´ files, normalize the counts using DESeq2.This requires a few steps:

    1. Install packages DESeq2, EdgeR
    1. Build a file description fileDesc.txt file for all .txt files about counts matrix,including file path, inoculation, and description.
    1. Read data files.
    1. Build meta data. meta file is similar to filedesc, bnut we use meta for create DESeqDataset object. Ensure the row names of the metadata dataframe are present and in the same order as the column names of the counts data frame.
    1. Create DESeqDataset object
    1. Generate size factors
    1. Generate the normalized counts

2.2.1 Install packages

Install DESeq2 package

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DESeq2")

We need to install EdgeR package because we need the readDGE function here.

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("edgeR")

2.2.2 Build a file description

Here the file description called fileDesc is already build and save in data folder. We will use it for reading a list of seperated files.

2.2.3 Read data files

Set to the current work directory

setwd('C:/Users/timoi/OneDrive - University of Helsinki/ZilanWen Github/RNA-seq_Differential-Gene-Expression-Analysis')

Now we start to read raw data

Coinfection.targets<-read.delim("./data/fileDesc.txt")

change the rawnames of the dataframe Coinfection.targets, later on we don’t need to rename column names.

rownames(Coinfection.targets)<-c("Ha1","Ha2","Ha3","Ctr1","Ctr2","Ctr3")

load the packege edgeR

library(edgeR)
## Loading required package: limma

read the six .txt files listed in Coinfection.targets.

Coinfection.orig <- readDGE(Coinfection.targets, header=F)

check the dimension of the data set

dim(Coinfection.orig)
## [1] 51751     6

chreck the first 6 rows of the data

head(Coinfection.orig)
## An object of class "DGEList"
## $samples
##           files.path inoculation         description           files group
## Ha1   ./data/Ha1.txt          Ha H.annosum infection  ./data/Ha1.txt     1
## Ha2   ./data/Ha2.txt          Ha H.annosum infection  ./data/Ha2.txt     1
## Ha3   ./data/Ha3.txt          Ha H.annosum infection  ./data/Ha3.txt     1
## Ctr1 ./data/Ctr1.txt         Ctr             Control ./data/Ctr1.txt     1
## Ctr2 ./data/Ctr2.txt         Ctr             Control ./data/Ctr2.txt     1
## Ctr3 ./data/Ctr3.txt         Ctr             Control ./data/Ctr3.txt     1
##      lib.size norm.factors
## Ha1  12510293            1
## Ha2   8033676            1
## Ha3  10094630            1
## Ctr1 10763763            1
## Ctr2 10539029            1
## Ctr3 12251172            1
## 
## $counts
##             Samples
## Tags         Ha1 Ha2 Ha3 Ctr1 Ctr2 Ctr3
##   PITA_19277 274 169 195   90   80  329
##   PITA_36893   0   0   0    0    0    0
##   PITA_43071 157 101 150  128  133  127
##   PITA_31484   0   0   0    0    0    0
##   PITA_22913   0   0   0    0    0    0
##   PITA_43120  62  42  22   53   24   15

Extact counts dataframe

Coinfection.rawCount <- Coinfection.orig$count
dim(Coinfection.rawCount)
## [1] 51751     6
head(Coinfection.rawCount)
##             Samples
## Tags         Ha1 Ha2 Ha3 Ctr1 Ctr2 Ctr3
##   PITA_19277 274 169 195   90   80  329
##   PITA_36893   0   0   0    0    0    0
##   PITA_43071 157 101 150  128  133  127
##   PITA_31484   0   0   0    0    0    0
##   PITA_22913   0   0   0    0    0    0
##   PITA_43120  62  42  22   53   24   15

2.2.4 Build meta data

We define sampletype: we only have two types of samples, Control and Ha-infection

sampletype <- factor(c(rep("Ha",3), rep("Ctr", 3)))

Build meta data frame

meta <- data.frame(sampletype, row.names = colnames(Coinfection.orig$count))

check the column name of counts dataframe

colnames(Coinfection.orig$count)
## [1] "Ha1"  "Ha2"  "Ha3"  "Ctr1" "Ctr2" "Ctr3"

check the rowname of meta dataframe

rownames(meta)
## [1] "Ha1"  "Ha2"  "Ha3"  "Ctr1" "Ctr2" "Ctr3"

Check that sample names match in both files

all(colnames(Coinfection.orig$count) %in% rownames(meta))
## [1] TRUE

2.2.5 Create DESeqDataset object

load the package DESeq2 Create DESEq2 dataset dds

library(DESeq2)
dds <- DESeqDataSetFromMatrix(Coinfection.orig, colData = meta, design = ~ sampletype)

If you want to view dds, you can use the code:

head(counts(dds))
##             Samples
## Tags         Ha1 Ha2 Ha3 Ctr1 Ctr2 Ctr3
##   PITA_19277 274 169 195   90   80  329
##   PITA_36893   0   0   0    0    0    0
##   PITA_43071 157 101 150  128  133  127
##   PITA_31484   0   0   0    0    0    0
##   PITA_22913   0   0   0    0    0    0
##   PITA_43120  62  42  22   53   24   15

2.2.6 Generate size factors

To perform the median of ratios method of normalization, DESeq2 has a single estimateSizeFactors() function that will generate size factors.

dds <- estimateSizeFactors(dds)

take a look at the normalization factor applied to each sample

sizeFactors(dds)
##       Ha1       Ha2       Ha3      Ctr1      Ctr2      Ctr3 
## 1.1734464 0.7897579 0.9805609 1.0611942 0.9794690 1.1512821

2.2.7 Generate the normalized counts

Now, to retrieve the normalized counts matrix from dds, we use the counts() function and add the argument normalized=TRUE.

normalized_counts <- counts(dds, normalized=TRUE)

Save the normalized_countsinto results to your local path.

write.csv(normalized_counts, file="./results/coinfection_normalized_counts_DESeq2.csv")

3 Sample-level quality Control

3.1 Principal Component Analysis (PCA)

3.1.1 PCA overview**

A useful initial step in an RNA-seq analysis is often to assess overall similarity between samples:

  • Which samples are similar to each other, which are different?

  • Does this fit to the expectation from the experiment’s design?

  • What are the major sources of variation in the dataset?

To explore the similarity of our samples, we will be performing sample-level QC using Principal Component Analysis (PCA) and hierarchical clustering methods. These methods/tools allow us to check how well similar the replicates are to each other (clustering) and to make sure that the experimental condition is the major source of variation in the data. Sample-level QC can also help identify any samples behaving like outliers; we can further explore any potential outliers to determine whether they need to be removed prior to DE analysis.

These unsupervised clustering methods are run using log2 transformed normalized counts. The log2 transformation improves the distances/clustering for visualization. Instead of using an ordinary log2 transform, we will be using regularized log transform (rlog), to avoid any bias from the abundance of low-count genes; Note1 below explains this in more detail.

For more details, please go to DEG workshop Salmon

3.1.2 Using our dataset for PCA

  • Transform normalized counts for the dataset

To improve the distances/clustering for the PCA and heirarchical clustering visualization methods, we need to moderate the variance across the mean by applying the rlog transformation to the normalized counts.

Transform counts for data visualization

rld <- rlog(dds, blind=TRUE)

The blind=TRUE argument is to make sure that the rlog() function does not take our sample groups into account - i.e. does the transformation in an unbiased manner. When performing quality assessment, it is important to include this option. The DESeq2 vignette has more details about this.

The rlog() function returns a DESeqTransform object, another type of DESeq-specific object. The reason you don’t just get a matrix of transformed values is because all of the parameters (i.e. size factors) that went into computing the rlog transform are stored in that object. We use this object to plot the PCA and heirarchical clustering figures for quality assessment.

  • start with PCA

The function plotPCA() requires two arguments as input: a DESeqTransform object and the “intgroup” (interesting group), i.e. the name of the column in our metadata that has information about the experimental sample groups.

plotPCA(rld, intgroup="sampletype")
## using ntop=500 top features by variance

If we want to save PDF file on own computer:

pdf("./results/PlotPCA_dds.pdf")
plotPCA(rld, intgroup="sampletype")
dev.off()

  • Exercise:

    1. What does the above plot tell you about the similarity of samples?
    1. Does it fit the expectation from the experimental design?
    1. What do you think the %variance information (in the axes titles) tell you about the data in the context of the PCA?

3.2 Hierarchical Clustering Heatmap

Similar to PCA, hierarchical clustering is another, complementary, method for identifying strong patterns in a dataset and potential outliers. The heatmap displays the correlation of gene expression for all pairwise combinations of samples in the dataset. Since the majority of genes are not differentially expressed, samples generally have high correlations with each other (values higher than 0.80). Samples below 0.80 may indicate an outlier in your data and/or sample contamination.

There is no built-in function in DESeq2 for plotting the heatmap for diplaying the pairwise correlation between all the samples and the hierarchical clustering information; we will use the pheatmap() function from the pheatmap package. This function cannot use the DESeqTransform object as input, but requires a matrix or dataframe. So, the first thing to do is retrieve that information from the rld object using a function called assay() (from the SummarizedExperiment package) that converts the data in a DESeqTransform object to a simple 2-dimensional data structure (a matrix in this case).

Extract the rlog matrix from the object using assay function which is part of the SummarizedExperiment package which is a DESeq2 dependency and is loaded with the DESeq2 library

rld_mat <- assay(rld)

Next, we need to compute the pairwise correlation values for all the samples. We can do this using the cor() function:

rld_cor <- cor(rld_mat)

check the output of cor(), make note of the row names and column names

head(rld_cor)
##            Ha1       Ha2       Ha3      Ctr1      Ctr2      Ctr3
## Ha1  1.0000000 0.9974935 0.9978500 0.9956639 0.9956698 0.9937710
## Ha2  0.9974935 1.0000000 0.9978706 0.9961065 0.9957914 0.9953450
## Ha3  0.9978500 0.9978706 1.0000000 0.9956039 0.9954369 0.9943809
## Ctr1 0.9956639 0.9961065 0.9956039 1.0000000 0.9967852 0.9946562
## Ctr2 0.9956698 0.9957914 0.9954369 0.9967852 1.0000000 0.9943313
## Ctr3 0.9937710 0.9953450 0.9943809 0.9946562 0.9943313 1.0000000

We will use meta for annotation argument later, let’s check meta :

head(meta)
##      sampletype
## Ha1          Ha
## Ha2          Ha
## Ha3          Ha
## Ctr1        Ctr
## Ctr2        Ctr
## Ctr3        Ctr

Install pheatmap package and load the package

install.packages("pheatmap")

?pheatmap # help page for pheatmap

Plot heatmap using the correlation matrix and the metadata object

library(pheatmap)
pheatmap(rld_cor, annotation = meta)

If we want to change color, we can try the codes below: R color in Brewer’s palettes can be available the link https://r-graph-gallery.com/38-rcolorbrewers-palettes

pheatmap function instruction are also available [pheatmap] (https://www.rdocumentation.org/packages/pheatmap/versions/1.0.12/topics/pheatmap)

heat.colors <- RColorBrewer::brewer.pal(6, "Blues")
pheatmap(rld_cor, annotation = meta, color = heat.colors, border_color=NA, fontsize = 10, 
        fontsize_row = 10, height=20)

If you want to save the heatmap in local path in your computer:

pdf("./results/PlotHeatmap_dds.pdf")
heat.colors <- RColorBrewer::brewer.pal(6, "Blues")
pheatmap(rld_cor, annotation = meta, color = heat.colors, border_color=NA, fontsize = 10, 
        fontsize_row = 10, height=20)

4 Differential expression analysis (DEA) using EdgeR

Differential expression analysis can be carried out using DESeq2 or using EdgeR. The DESeq2 methods has been describes in DGE workshop Salmon. Here we will learn to use another method – EdgeR

Actual raw integer read counts (un-normalized) are then used for DGE analysis using edgeR. edgeR prefers the raw integer read counts.

To know the current directory

getwd()

Navigate to your work directory, for example: At RStudio –> Session –> Set Working Directory –> Choose Directory –> Desktop –> RNA-seq_DEA

or use code below:

setwd("C:/Users/username/Desktop/RNA-seq_DEA")

load EdgeR package

library(edgeR)
options(digits=3)

To tells R where the data files are

infection.targets<-read.delim("./data/fileDesc.txt")

To check Coinfection.targets

infection.targets
##        files.path inoculation         description
## 1  ./data/Ha1.txt          Ha H.annosum infection
## 2  ./data/Ha2.txt          Ha H.annosum infection
## 3  ./data/Ha3.txt          Ha H.annosum infection
## 4 ./data/Ctr1.txt         Ctr             Control
## 5 ./data/Ctr2.txt         Ctr             Control
## 6 ./data/Ctr3.txt         Ctr             Control

change the rawnames of the dataframe Coinfection.targets, later on we don’t need to rename column names.

rownames(infection.targets)<-c("Ha1","Ha2","Ha3","Ctr1","Ctr2","Ctr3")

To check Coinfection.targets again, The difference with the above one is an addtional rownames on the left side

infection.targets
##           files.path inoculation         description
## Ha1   ./data/Ha1.txt          Ha H.annosum infection
## Ha2   ./data/Ha2.txt          Ha H.annosum infection
## Ha3   ./data/Ha3.txt          Ha H.annosum infection
## Ctr1 ./data/Ctr1.txt         Ctr             Control
## Ctr2 ./data/Ctr2.txt         Ctr             Control
## Ctr3 ./data/Ctr3.txt         Ctr             Control

read and merges a set of text files containing gene expression counts, it makes a DGEList object directly

infection <- readDGE(Coinfection.targets, header=F)

Check the dimension of DGElist R object

dim(infection)
## [1] 51751     6
head(infection)
## An object of class "DGEList"
## $samples
##           files.path inoculation         description           files group
## Ha1   ./data/Ha1.txt          Ha H.annosum infection  ./data/Ha1.txt     1
## Ha2   ./data/Ha2.txt          Ha H.annosum infection  ./data/Ha2.txt     1
## Ha3   ./data/Ha3.txt          Ha H.annosum infection  ./data/Ha3.txt     1
## Ctr1 ./data/Ctr1.txt         Ctr             Control ./data/Ctr1.txt     1
## Ctr2 ./data/Ctr2.txt         Ctr             Control ./data/Ctr2.txt     1
## Ctr3 ./data/Ctr3.txt         Ctr             Control ./data/Ctr3.txt     1
##      lib.size norm.factors
## Ha1  12510293            1
## Ha2   8033676            1
## Ha3  10094630            1
## Ctr1 10763763            1
## Ctr2 10539029            1
## Ctr3 12251172            1
## 
## $counts
##             Samples
## Tags         Ha1 Ha2 Ha3 Ctr1 Ctr2 Ctr3
##   PITA_19277 274 169 195   90   80  329
##   PITA_36893   0   0   0    0    0    0
##   PITA_43071 157 101 150  128  133  127
##   PITA_31484   0   0   0    0    0    0
##   PITA_22913   0   0   0    0    0    0
##   PITA_43120  62  42  22   53   24   15

Get the raw mapped count before filtering

infection.rawCount <- infection$count

Export raw count table into results folder

write.csv(infection.rawCount, file="./results/infection.rawCounts.csv")

Get the counts per million (TMM normalised) before filtering

infection.normCPM <- cpm(calcNormFactors(infection))
dim(infection.normCPM)
## [1] 51751     6
head(infection.normCPM)
##             Samples
## Tags           Ha1   Ha2   Ha3  Ctr1  Ctr2  Ctr3
##   PITA_19277 22.35 20.46 19.04  8.12  7.85 27.31
##   PITA_36893  0.00  0.00  0.00  0.00  0.00  0.00
##   PITA_43071 12.81 12.23 14.65 11.55 13.06 10.54
##   PITA_31484  0.00  0.00  0.00  0.00  0.00  0.00
##   PITA_22913  0.00  0.00  0.00  0.00  0.00  0.00
##   PITA_43120  5.06  5.09  2.15  4.78  2.36  1.24
write.csv(infection.normCPM, file="./results/infection.normCPM.csv")

Keep genes that are expressed at least 1 CPM in at least 3 libraries, normally it is the number of biological replicates of smaller group

infection.filtered <- rowSums(cpm(infection)>1) >=3
table(infection.filtered)
## infection.filtered
## FALSE  TRUE 
## 25742 26009

Libraries size of data BEFORE filtering

infection$samples$lib.size
## [1] 12510293  8033676 10094630 10763763 10539029 12251172

cover the original file with our filter data

Infection <- infection[infection.filtered,]

libraries size of data after filtering

colSums(Infection$counts)
##      Ha1      Ha2      Ha3     Ctr1     Ctr2     Ctr3 
## 12455638  7998136 10047983 10717381 10497878 12206822
dim(Infection)
## [1] 26009     6

Update the filtered libraries size

Infection$samples$lib.size <- colSums(Infection$counts)
Infection$samples
##           files.path inoculation         description           files group
## Ha1   ./data/Ha1.txt          Ha H.annosum infection  ./data/Ha1.txt     1
## Ha2   ./data/Ha2.txt          Ha H.annosum infection  ./data/Ha2.txt     1
## Ha3   ./data/Ha3.txt          Ha H.annosum infection  ./data/Ha3.txt     1
## Ctr1 ./data/Ctr1.txt         Ctr             Control ./data/Ctr1.txt     1
## Ctr2 ./data/Ctr2.txt         Ctr             Control ./data/Ctr2.txt     1
## Ctr3 ./data/Ctr3.txt         Ctr             Control ./data/Ctr3.txt     1
##      lib.size norm.factors
## Ha1  12455638            1
## Ha2   7998136            1
## Ha3  10047983            1
## Ctr1 10717381            1
## Ctr2 10497878            1
## Ctr3 12206822            1

Performed normalisation with TMM method

Infection = calcNormFactors(Infection)

The libraries after normalisation

Infection$samples
##           files.path inoculation         description           files group
## Ha1   ./data/Ha1.txt          Ha H.annosum infection  ./data/Ha1.txt     1
## Ha2   ./data/Ha2.txt          Ha H.annosum infection  ./data/Ha2.txt     1
## Ha3   ./data/Ha3.txt          Ha H.annosum infection  ./data/Ha3.txt     1
## Ctr1 ./data/Ctr1.txt         Ctr             Control ./data/Ctr1.txt     1
## Ctr2 ./data/Ctr2.txt         Ctr             Control ./data/Ctr2.txt     1
## Ctr3 ./data/Ctr3.txt         Ctr             Control ./data/Ctr3.txt     1
##      lib.size norm.factors
## Ha1  12455638        0.981
## Ha2   7998136        1.021
## Ha3  10047983        1.011
## Ctr1 10717381        1.031
## Ctr2 10497878        0.969
## Ctr3 12206822        0.988

Get the counts per million (TMM normalised) after filtering

Infection.filtered.normCPM <-cpm(calcNormFactors(Infection))

Export TMM normalized count table after filtering

write.csv(Infection.filtered.normCPM, file="./results/Infection.filtered.normCPM.csv")

differentailly gene expression analysis

Treatment factor

group<-factor(c('Ha','Ha','Ha',"Ctr","Ctr","Ctr"))

Describe the experimental design,one factor with intercept, here Ctr is the intercept

Infection.design <- model.matrix(~group)   
rownames(Infection.design)<-colnames(Infection$counts)
Infection.design
##      (Intercept) groupHa
## Ha1            1       1
## Ha2            1       1
## Ha3            1       1
## Ctr1           1       0
## Ctr2           1       0
## Ctr3           1       0
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"

to see if some samples are outliner

plotMDS(Infection, main="MDS plot of RNA-Seq", labels=colnames(Infection$counts))

Export the MDS plot to a PNG file

png("./results/plotMDS.Infection.png", res=300, height=1800, width=1800)
plotMDS(Infection, main="MDS plot of Infection RNA-Seq", labels=colnames(Infection$counts))
dev.off()

Estimating the common dispersion, trended dispersion, tagwwise dispersion Dispersion means biological coeffient of variation (BCV) squared. E.g. if genes expression typically differs from replicate to replicate by 20% its BCV is 0.2, and its dispersion is 0.04. Estimating the common dispersion

Infection <- estimateGLMCommonDisp(Infection, Infection.design)

Estimating the trended dispersion

Infection <- estimateGLMTrendedDisp(Infection, Infection.design)

Estimating the tagwwise dispersion

Infection <- estimateGLMTagwiseDisp(Infection, Infection.design)

Creating a visual representation of the mean-variance relationship and BCV-mean logCPM relationship

plotMeanVar(Infection, show.tagwise.vars=T,NBline=T)
plotBCV(Infection)